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Abstract 

We have investigated a single charged microgel in aqueous solution with a combined simulational 
model and Poisson-Boltzmann theory. In the simulations we use a coarse-grained charged bead- 
spring model in a dielectric continuum, with explicit counterions and full electrostatic interactions 
under periodic and non-periodic boundary conditions. The Poisson-Boltzmann model is that of a 
single charged colloid confined to a spherical cell where the counterions are allowed to enter the 
uniformly charged sphere. We compare the simulational results to those of the Poisson-Boltzmann 
solution and find good agreement, i.e., for the number of confined counterions within the gel. We 
then proceed to investigate the origin of the differences between the results these two models give, 
and performed a variety of simulations which were designed to test for the influence of charge 
correlations, excluded volume interactions, and thermal fluctuations in the strands of the gel. Our 
results support the applicability of the Poisson-Boltzman cell model to study ionic properties of 
small microgels under dilute conditions. 
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I. INTRODUCTION 



Hydrogel research has constantly attracted a lot of attention, due to its capability of ab- 
sorbing large quantities of water, resulting in the expansion of the material up to 1000 
times its dry volume. Recent review articles on hydrogel research range from general 
reviews- to application specific reviews such as sensors^, surface applications^., drug delivery 
applications^ and other biomedical applications^i^i^. 

A hydrogel is a crosslinked network of water soluble polymers. The hydrophilicity of 
the free uncrosslinked chain allows it to attract water molecules and be soluble in water up 
to infinite dilution. The crosslinking, however, restricts the dilution to a maximal swelling 
of the gel. The combination of these two properties allows the hydrogel to attract water 
molecules and keep them within the network, resulting in an overall swelling of the entire 
gel. Hydrogels can be made of neutral or charged polymers called polyelectrolytes (the 
backbone of polyelectrolytes, however, need not be hydrophilic) . When polyelectrolyte gels 
are dissolved in water, charges are formed in the polyelectrolyte when the counter ions 
dissociate from the polymer. These counter ions are free to move within the solvent, but 
stay, in the absence of salt, within the vicinity of the polyelectrolyte to keep electrostatic 
neutrality of the gel. The increase in entropy due to these counter ions overcomes the extra 
enthalpy to bring in even more water molecules into the network. The overall effect is that 
charged hydrogels have an even greater swelling ratio than their uncharged counterparts. A 
hydrogel matrix can be made to a few centimeters in diameter, depending on the synthetic 
procedure used. Microgels, on the other hand, are hydrogels that are within the micrometer 
range. 

A number of analytical theories on hydrogels and polyelectrolytes have been 
presented^'^i^'iiii^. Corresponding simulations have also been done^^^^^, confirming the 
basic picture that the entropy of the counter ions exerts an extra osmotic pressure which 
in turn increases the swelling capacity of the hydrogels. Thus determining the location of 
these counter ions within the network, and specifically the investigation of ion condensation 
around the polymer chains, is an important aspect of hydrogel research. To our knowledge, 
the simulations that have been done, however, all consider quasi infinite systems via periodic 
boundary conditions that connect the chains of one side of the simulation box to its opposite 
side. This set up suffers from the limitation that the counter ions are, in a sense, always 
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inside the polymer matrix. Thus they cannot examine the actual charge fraction of counter 
ions inside the network which could, in turn, change the amount of osmotic pressure exerted 
by the counter ions, and thus change the swelling properties. 

In an attempt to attack this problem, we simulate a model microgel, which is a finite 
sized object made up of crosslinked polymers in a periodic box. The simulation details are 
described in section III A[ The system is a spherical polymer network in a simulation box, 
as shown in figure [TJ 

This model simulation allows us to calculate the counter ion condensation not only within 
the vicinity of each chain, but also the fraction inside the sphere, and study the consequences 
for the chain conformations. This model is also similar to an even more reduced model 
of a single charged sphere with counter ions that can be situated inside or oustide this 
charged sphere. This reduced model can be calculated using Poisson-Boltzmann (PB) theory. 
The simulation results can therefore be compared to a PB calculation, and thus provide 
evidence for the applicability of the Poisson-Boltzman cell model to study ionic properties 
of microgels. 

There exists, however, quite a number of differences between the Poisson-Boltzmann 
method and simulation set-up, which can potentially confuse the interpretation of a direct 
comparison between the two. This paper examines the effects of these differences by pre- 
senting model simulations aimed at systematically removing, or adding, these factors in 
a step-wise fashion, with the final aim of studying the usefulness of the PB theory for a 
description of the swelling properties of microgels. 

In section ini we describe the details of the microgel (MG) simulation (III Al) . brief descrip- 
tion of PB theory (III BP . and lastly the differences between PB calculation and the microgel 
simulation fill CI) . This list of differences (or assumptions) provides the rationale for the rest 
of the model simulations that we performed, aimed at connecting, in a stepwise manner, 
the PB model to the MG simulation. The methodologies and results of these simulations 
described in sections [Tlll through |Vl And in section IVll we give our counclusions. 
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II. MICROGEL SIMULATION AND POISSON-BOLTZMANN CALCULATION 
METHODOLOGIES 



In all simulations, the following reduced units were used: energy is placed in terms of 
ksT which is given a value of 1.0, values of lengths are in a, charge is in units of e, and time 
is in r. All particles are given a mass of 1.0. Table [T] summarizes the parameters that are 
the same in all simulations in this paper, whereas table [III and IIIII lists the parameters that 
are specific to each simulation. All simulations were done using ESPResSoiS., a molecular 
dynamics package for soft matter simulations, under NVT conditions. In order to simulate 
a constant temperature ensemble, a Langevin thermostat is used, where all particles are 
coupled to a heat bath. Thus the equation of motion is modified to 

= -V^Vtot (r1) - m,r^ + (t) (1) 

where rrii is the mass of particle i, Vtot is the total potential acting on the particle, F is the 
friction coefficient that is set to 1.0, and Wi{t) is a stochastic force that is uncorrelated in 
time and across particles. The friction and stochastic forces are linked by the fluctuation- 
dissipation theorem 

(t) ■ W, (t)) = QmYksTd^jd {t - t') (2) 

where = 1.0. 



A. Microgel Simulation Details 

The parameters for the microgel simulation are listed in Tables H] and llllt labelled as 
P-MG, where P indicates that periodic boundary conditions were used. Figure [T^ shows 
how the starting structure of the microgel was constructed, whereas figure gives a snap- 
shot of an equilibrated microgel configuration. The microgel was simulated using a general 
bead-spring modeir^^. The microgel was composed of 46 polymers, each composed of 10 
monomers (beads), giving a total of 460 monomers. These polymers were connected by 
29 crosslinks, each represented by one bead, giving a total of 489 beads (monomers plus 
crosslinks) in the microgel. Some crosslinks, especially those inside the microgel, were con- 
nected to the ends of four different polymers (serving as tetrafunctional nodes), while the 
rest were conneced to the ends of three different polymers (trifunctional nodes). Since all 
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of the polymer ends were connected to one crosslink, the microgel had no dangling chains. 
Upon examining figure [1^, one could see that the starting structure resembles that of a 
diamond lattice, with the edges cut off to approach a spherical shape. Once the particles are 
allowed to move, the average shape further approximates that of a sphere. The whole set-up 
was chosen to have an average equilibrated shperical shape, in order to make a comparison 
with the PB calculation of colloid using a cell model. 

Half of the monomers (that is, 230) were given a charge of +1.0, while the others remained 
uncharged. The charged and uncharged monomers were placed along the polymer in an 
alternating sequence. Out of the 29 total crosslinks, 20 were also given a charge of +1.0, 
resulting to a total of 250 of the 489 total particles in the microgel having a charge of +1.0. 
Thus 250 free counter ions, each with a charge of -1.0, were added to the simulation, giving 
a total 739 objects and an overall neutral charge. 

The non-bonded interactions of the particles were described by a Weeks- Chandler- 
Anderson potential^ of the form 

I 4er., ( ^1 

ULjirij) 
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(3) 

for Tij > Tcut 



where (Tq = 1.0, eij = 1.0, Vcut = 2^^^, and Cghift = 0.25, thus leaving only the repulsive part 
of the interaction. These are the standard parameters for bead-spring model of polymers in 
a good solvent^^. The bonds between the particles were modelled using the FENE (Finite 
Extension Nonlinear Elastic) potential^ 
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for r < rp 

(4) 

oo for r > rp 



where kp = 10.0^^ and rp = 1.5^. These parameters allow the bonds to equlibrate to 
r ^ 1.0. 

The potential energy between charges and qj separated by a distance r^j were described 
by the unscreened Coulomb potential 

Ucinj) = — (5) 

I ij 

where Ib, the Bjerrum length, is defined as 



where cq is the elementary charge, eo is the vacuum permittivity, and e,. is the relative 
dielectric constant of the solvent. The Bjerrum length is defined as the distance at which 
the Coulomb energy between two unit charges is equal to the interaction energy of ksT. This 
parameter provides a convenient way of tuning the strength of electrostatic interactions in 
a simulation. As a point of reference, the Bjerrum length of water at room temperature is 
7.1 A. In this paper, a Bjerrum length of /b = 1 and ksT = 1 was used in all simulations 
(Table [T]). Given a typical hydrogel such as poly(N-isopropylacrylamide) (poly-NIPAAm) , 
the average distance between the centers of masses between two repeat units is within the 
order of 7.1 A, making the choice of r = 1.0 and Ib = 1-0 both computational convenient 
and comparable to values of actual hydrogels. 

The electrostatics for this simulation box was treated using the P3M algorithmSfi'H, tuned 
to an accuracy of 10~^ in the absolute error of the Coulomb forces with the help of the 
formulas in Re^. 

The MG simulation was performed in a cubic box with a length of 84 and periodic 
boundary conditions. This size was determined to be much larger than the Debye length of 
the system. The Debye length is the screening distance of a charged particle in a solution, 
beyond which its interaction with other ions is effectively screened. The Debye length is 
given as 

where hq is the molar conentration of ions in the reservoir, that is, the concentration of 
the counter ions remaining outside the sphere. As will be seen in the results, the microgel 
equilibrates to a radius of 15, with an average of 46% of the counter ions remaining outside 
the micorgel. This results to a Debye length of 14.1. The closest distance between two 
periodic images of the microgel is 54, which is roughly 3.8 times this Debye length. This 
method of calculating the molar concentration uq already yeilds the upper limit of the 
Debyle length. Choosing any other method for defining the Uq — for example, Uq as the total 
number of counter ions over the total box volume — would yeild a higher value for uq, and 
thus a shorter Debye length. This set-up therefore guarantees enough screening to prevent 
the microgel from being strongly electrostatically affected by its periodic image. Thus the 
microgel is in the dilute limit where each micgrogel behaves independently of other microgels. 
This allows us to compare the results with the PB calculations of a colloid in a cell model, 
which also assumes that the interaction of one colloid with other colloids is negligible. 
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B. Iterative Poisson-Boltzmann Solution 

The Poisson-Boltzmann (PB) theory is a mean-field theory that describes the electrostatic 
interactions between charged particles in solution. The Poisson-Boltzmann equation is 



where 0(r) represents the electrostatic potential, and Zi is the valency of the ion species i. 
The charge density p(r), assumed to be a Bolztmann distribution of charges, is given by 



This theory assumes that the only interactions taken into consideration are Coulombic in- 
teractions between point charges, thus neglecting any finite size effects of the charges. The 
other assumptions are that the ion entropy is that of an ideal gas. 

In a cell model^^, one can approximate the interactions of a charged object with its counter 
ions in a local volume, and assume that the inteaction of this unit with other similar units is 
negligible. Using this model, we can therefore think of a microgel as one such charged object 
isolated in a spherical volume. This spherical cell model is shown in in figure [2^. The inner 
circle represents the charged object — for example, a colloid or a microgel — where the total 
(positive) charge of the object is distributed evenly within the sphere. This charged object 
has radius tq. It is co centric to an outer spherical cell of radius R, as represented by the 
outer circle in the figure. Counter ions are present in this spherical cell. In this study, these 
counter ions were represented as point charges with a charge of —1.0. They were allowed to 
go anywhere within the outer sphere, even into the charged object. 

Our MG set-up (figures [1] and [2^) is fairly close to a charged sphere inside a spherical 
cell model (figure [2|i). Solving the PB equation for this model of a charged sphere would 
therefore be a good comparison for the resulting distribution of charges of the MG simulation. 
The PB equation is a non-linear partial differential equation, and has closed-form analytical 
solutions only for a limited number of cases. For our case, it has to be solved numerically, 
as was done by Deserno et. al.—'^, where the PB equation was solved numerically for the 
same spherical cell model as described above. 

Using this PB numerical solver, we calculated the fraction of counter ions present in 
the charged sphere. The radius of the outer spherical cell was set to i? = 52. The inner 
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charged sphere had a radius of tq = 18 and was given a total charge of +250.0. There 
were nci = 250 counter ions present, each with a charge of —1.0. The Bjerrum length was 
chosen to be Ib = 1-0. These parameters were chosen so as to match those of the MG 
simulation: the volume of the spherical cell is equal to the volume of the simulation box, 
and the charged sphere radius of rg = 18 is equal to the energy minimized structure (figure 
^) of the microgel coming from the initial structure (figure [1^). The parameters and results 
of this calculation, labelled as PB-1, are listed in Table HTl The results will be discussed in 
the suceeding sections. One other PB calculation was done, labelled as PB-2, with slightly 
different parameters listed in IIIII The reason for this will be explained in the discussion of 
the MG simulation in section |Vl 

C. Rationale of the Simulation Strategy 

The charged sphere in a cell model is already quite close to the microgel simulation. 
However, there still are quite a number of assumptions made between the PB model and 
the MG simulation. These are outlined as follows: 

1. Discrete Particles. Whereas the total charge of the sphere in the PB model is 
uniformly distributed throughout the charged sphere, discrete charges are used in the 
MG simulation. 

2. Excluded Volume. The microgel particles, whether charged or uncharged, in the MG 
simulation have excluded volume, as parameterized by ctlj, which in turn dertermines 
the local binding energy. This is non-existent in the charged sphere of the PB solver. 

3. Arranged Particles. The charged particles in the MG simulation are situated along 
the chain, and are therefore not totally free to move relative to each other, differing 
from the even distribution of charges in the PB model. 

4. Thermal Fluctuations. Aside from being arranged along the chains in the microgel, 
the charged particles also move as the polymer itself moves due to thermal fluctuations. 
This is not included in the PB theory. 

5. Periodic Boundary Conditions. Lastly, the PB model was calculated in an isolated 
cell model, whereas the MG simulations were done under periodic boundary conditions. 
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This difference also compares the electrostatic algorithms employed in the two systems. 
When a simulation is done in a cell model, a direct sum of the Coulombic potential 
between all pairs of charged particles is calculated, whereas an Ewald sum approach — 
in this paper, the P3M algorithm — is done for the MG simulation in the periodic 
box. 

Due to these assumptions, one would easily suspect that there would also be differences 
between the results of the PB calculation and the microgel simulation. However, it would 
be impossible to directly link the differences in results to each of these assumptions, since 
all of these assumptions are simultaneously present and could easily confuse the analysis. 
We therefore performed a series of model simulations that systematically remove from the 
original model one or more of the assumptions listed, thereby going from the PB picture to 
the MG simulation in a step-wise fashion. Using this procedure, we are able to attribute 
specific effects to the corresponding assumptions. These simulations are called Lattice Sphere 
(C-LS), Static Microgel (C-ST), and Microgel in a Sperical Shell (C-MG), which will be 
discussed below. The C in the labels indicate that they are done using a cell model, as 
compared to the microgel (P-MG), which was simulated using periodic boundary conditions 
(hence, labelled P). The set-ups of these simulations are shown in figure [21 The numbers on 
the arrows in this figure indicate the assumptions removed when going from one simulation 
to the next. 

In all the calculations and simulations in this paper (as shown in figure [2]), the following 
variables were kept constant: 

1. Total volume. The total volume in which counter ions where free to move were 
the same for all simulations and PB calculations. An outer radius of R = 52.0 was 
used for all the calculations involving the spherical cell model. A cube with the same 
volume (within a 0.63% difference) has a length of L = 84, which is the value used for 
the simulations invovling a cubic box. 

2. Counter Ions. The total number of counter ions was rici = 250, each with a charge 
of qci = -1.0. 

3. Chcirge of Inner Sphere. In order to keep the system electrically neutral, the total 
charge of the inner sphere was kept at g^p = +250.0. In some simulations, however, 
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the total number of charged particles in the inner sphere were greater than 250. In 
these cases, the fractional charge of these particles was adjusted in order to keep the 
total charge at +250.0. 

III. LATTICE SPHERE SIMULATIONS (C-LS) 
A. C-LS Methodology 

We simulated a series of systems which we call Lattice Sphere (C-LS), which consist of 
charged particles that are fixed on a cubic grid. The set-up is shown in figure [2]d. The final 
shape of this lattice approximates a sphere. This spherical lattice of radius tq = 18 was 
fixed at the center of an outer sphere with radius R = 52. Counter ions were placed inside 
the outer sphere. The outer sphere is given the same LJ repulsive term with parameters 
in Table [T] so as to keep the counter ions inside the sphere. During the simulation, the 
counter ions are free to move throughout the whole sphere (note that i? = 52 is the radius 
of the volume available for counter ion motion). A repulsive Lennard- Jones potential was 
used between counter ions and the charged particles in the lattice so as to avoid that two 
oppositely charged particles to interact too strongly via Coulombic forces. No Lennard- 
Jones potential was used between counter ions, since these were treated as point particles, 
following the assumption of PB theory (section III B[) . This set-up is very similar to the 
PB model, except that two assumptions were already were lifted: discrete charged particles 
(assumption [1]) distributed as evenly as possible within a sphere, and the added excluded 
volume (assumption [2]) of these particles instead of an even charge distribution. 

All C-LS simulations had a total of rici = 250 counter ions each with a charge of —1.0, 
and a total charge of +250.0 for the lattice. Given the various possible cubic lattice arrange- 
ments of charged particles with an overall spherical shape, the three arrangements used 
in this study had 250, 484 and 894 charged particles in the lattice sphere. And, in order 
for this sphere to retain a total of +250, charges of +1.0, +0.516, and +0.279, were dis- 
tributed in these respective arrangments. In doing so, the point charges were systematically 
smeared out within the sphere, going from 250 (+1.0) to 484 (+0.516) and lastly to 894 
(+0.279), therefore lessening the effect of charge discretization. In addition, the ctlj was 
also decreased from a = 1.0 to smaller values. This has the effect of minimizing the effect 
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of excluded volume (assumption [2]). 

We used two ways of determining these smaller aij values. The first way was to maintain 
the energy of closest approach, that is, the Coulombic energy between the counter ion (with 
a charge of -1.0) and the charged particle (whether charged +1.0, +0.516, or +0.279) in 
the three different systems were always be the same. By examining equation O varying 
both the charge of the particle (g^ in the numerator) and the minimum distance (r^ in 
the denominator) in the same way would yield equal energies. Thus, when varying the the 
charge from +1.0 to +0.516 to +0.279, we also varied the a of the LJ potential in equation[3] 
from 0" = 1.0 to cr = 0.516 to a = 0.279. These are done in three calculations: C-LS-1 with 
a = 1.0, qi = +1.0, and ucp = 250; C-LS-2 with a = 0.516, = +0.516 and ucp = 484; 
and C-LS-3 with a = 0.279, = +0.279 and ricp = 894. Therefore the differences in 
results of these three calculations cannot be due to the energy of closest approach. 

The second procedure used in decreasing the values of cr is by maintaining the ratio of 
the volume occupied by the charged particle to the volume of the entire inner sphere. In 
the previous procedure, the ratio of volumes of the charge particle to the inner sphere were 
1.68 X 10^3 for C-LS-1, 4.44 x 10"^ for C-LS-2, and 1.30 x 10"^ for C-LS-3. Keeping 
this ratio equal to C-LS-1, the a^j was changed to the following values for a new set of 
calculations: C-LS-4 with au = 0.802, q^ = +0.516 and ncp = 484; and C-LS-5 with 
aij = 0.654, qi = +0.279, and ricp = 894. Thus having kept the free volume inside the 
inner sphere constant, we can safely assume that the differences in results of these three 
calculations would not be due to entropic effects, which have been significantly reduced in 
this procedure. 

These five Lattice Sphere simulations aimed at isolating the effects of the first two PB 
assumptions. The parameters and results are summarized in table [Tll These simulations 
were all done in a spherical cell model, with the electrostatics potential computed as a direct 
sum of Coulombic energies between all pairs of particles. 

B. Discussion: Result of Charge Discretization and Excluded Volume 

Refering to table [Tll the first discussion compares the results of PB-1, C-LS-1, C-LS-2, 
and C-LS-3. These four calculations have the following parameters that are the same: total 
charge of sphere, number of counter ions, radius of charged sphere/lattice, and outer radius. 
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From calculations C-LS-1 to C-LS-2 to C-LS-3, the total number of discrete charged 
particles in the inner sphere was increased from 250 to 484 to 894, with a corresponding 
decrease in partial charge and cxij from 1.0 to 0.517 to 0.279. Figure [3] shows the integrated 
counter ion probability of the simulations listed in Table [TTl The second to the last line in 
Table [III refer to the results of the percentage counter ion in figure [3] where r is equal to the 
sphere radius tq. The dots in the inset of figure [3] at rg = 18a indicate these percentages. 
Figure H] shows the density of the counter ions as a function of the distance from the center 
of the inner sphere. 

In general, the integrated counter ion probability (Fig. [3]) and the counter ion density 
(Fig. H]) graphs of the C-LS simulations are close to the PB-1 graph, with C-LS-1 being 
the farthest and lower and PB-1, C-LS-2 being closer but higher than PB-1, and C-LS-3 
being practically the same as the PB-1. The fraction of counter ions inside the sphere at 
tq = 18 (last line of Table [TTl and the inset in Figure [3]) also reflects this trend. As the 
number of charged particles is increased (with the corresponding decrease in partial charge) 
and thus the model becoming closer to the PB picture, the result of counter ion distribution 
approaches the PB calculation. Of the three LS calculations, the one with a unit charge of 
+ 1.0 can be considered closest of real physical systems, since charged particles always have 
unit charges. Thus we can say that the discretization of charges in a real system would have 
a counter ion condensation around 3% lower than what is predicted by PB theory. 

However, there seems to be a slight anomaly in the trend shown above. Although a 
stepwise smearing out of the charges should have corresponded to a gradual approach to 
the PB-1 calculation (52.83%), the results showed a jump from the lower value for C-LS-1 
(49.80%), to higher for C-LS-2 (53.85%), then back down to practically equal for C-LS-3 
(52.78%) to the PB-1 calculation. This trend cannot be caused by the energy of closest 
approach, since this has been set-up to be equal in these three LS calculations. We then 
proceed to examining the second assumption applied to the LS calculations, that of the 
excluded volume of the charged particles. The three calculations whose exculded volume, as 
measured by the ratio of volumes of charge particle to that of the inner sphere, are C-LS-1, 
C-LS-4, and C-LS-5. The integrated counter ion distribution and the counter ion densities 
of C-LS-2 and C-LS-4 are practically the same. The same is true for the graphs of C-LS-3 
and C-LS-5. The percentage of counter ions inside the inner sphere are listed at the bottom 
of table [TTl The trends are the same as the first set of LS calculations: C-LS-1 is lower 
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(49.80%) than PB-1 (52.83%), C-LS-4 is higher (53.35%), and C-LS-5 (52.23%) is closest 
to PB-1. This shows that the counter ion distributions is not dependent on the excluded 
volume of the charged particles. 

The explanation of this trend simply comes from the ground state energies of the different 
set-ups, since the ground-state structures depend mainly on the configuration of the system. 
To do this, we set the Langevin temperature T = 0.0 and using a small friction coefficient 
r between 0.01 to 0.001 to allow for slow relaxation to the ground state. This resulting 
trend of the ground state electrostatic energies of these systems are similar to the trends 
in percentage counter ion in the sphere: C-LS-1 being the highest —2160.72^57" (least 
stable), C-LS-2 being the highest — 2301.29/^5^ (most stable), and C-LS-3 in the middle 
— 2293.17fc5T. Thus the trend in counter ion fraction in the lattice sphere depends mostly 
on ground state configurations rather than the excluded volume of the charged particles. 
This does not mean, however, that excluded volume will not have an effect in polymeric 
systems. The next section discusses the more realistic arrangement of charged particles 
within a polymer chain. 

IV. STATIC MICROGEL SIMULATION DETAILS (C-ST) 
A. C-ST Methodology 

After having seen the effects of assumptions [1] and [2], we then proceeded to simulate 
another system with assumption [3] removed, that of charged particles in the system being 
placed along the chains in the microgel. In this case, we took the coordinates of an expanded 
microgel with the same parameters as described in III Al and fixed their coordinates so as 
not to allow thermal fluctuations. A snapshot of this system is shown in Figure Wp- This 
simulation is labelled as C-ST. Two versions of this simulation were performed. In the first 
case (C-ST-1), we removed the uncharged monomers (total of 230) so as not to include 
the effects of their excluded volume in the simulation. In the second case (C-ST-2), the 
uncharged monomers were not removed, so that the microgel has the same number of charged 
and uncharged particles present in the microgel simulations (C-MG and P-MG). These 
two therefore tested the effect of the arranging the charged particles along a chain, with the 
second looking at the added effect of the excluded volume of the uncharged particles, which 
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is therefore one step closer to the MG set-up. The parameters are similar to PB and C-LS, 
as listed in Table [Tll This particular arrangement of the particles had an overall effect of 
lessening the space between the charged particles beside each other on the same chain, and 
increasing empty spaces between charged particles on different chains. As with the Lattice 
Sphere case, this simulation was also done in a spherical cell model, with the electrostatics 
calculated direct sum. 

B. Discussion: Arrangement of Excluded Volume along a Polymer 

The parameters and results of the static microgel simulations, C-ST-1 and C-ST-2, are 
listed in Table HTj with the integrated counter ion probability curve shown in figure [3] and the 
counter ion density shown in figure |H The results show a further decrease of the fraction of 
counter ions in the static microgel C-ST-1 (46.28%) and in C-ST-2 (45.67%) as compared 
to PB-1 (52.83%) and C-LS-1 (49.80%). Figure [5] shows the relative distances of particles 
for the calculations of: a) C-LS-1, b) C-ST-1, and c) C-ST-2. The black circles represent 
the charged particles in the lattice sphere or microgel, the grey circle in fig. [5t represents 
the uncharged particle in C-ST-2, the dotted circles represent the excluded volume of these 
particles (Jlj-, and the unshaded circles represent the counter ions. 

In the case of the C-LS-1 configurations (fig. [5^), there is enough space in between 
two charged particles for a counter ions to move around the charged particle. The closest 
distance between the two dotted circles (given by x) is greater than 1.0, which is the Bjerrum 
length in the simulation. Thus even in this case, their Coulombic repulsion would still be 
much less than their thermal energy. 

That is not the case for the static microgel (C-ST-1, fig. [Sb) which does not include 
the uncharged particles in the microgel. The closest distance between paths (dotted circles) 
of two counter ions of adjacent charged particles is less than the Bjerrum length. Thus, 
unlike the lattice sphere cases, the space right in between two charged particles cannot be 
occupied by two counter ions, thus lessening the overall free space available for the counter 
ions. This effect could in turn decrease the fraction of charges in the microgel. The results of 
the second static microgel simulation (C-ST-2, fig. [5b) further supports this point. In this 
case, the additional uncharged particle (grey circle) contributes the extra excluded volume to 
the same space where electrostatic interactions prevented the counter ions to go, as argued 
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in the case of C-ST-1. Unlike the 3.52% drop from C-LS-1 to C-ST-1, these two static 
microgel simulations have very close percentages of counter ions inside the gel — 46.28% 
for C-ST-1 and 45.67% for C-ST-2 — indicating that the additional uncharged particles 
contributed only slightly to the reduction of counter ions in the gel. Thus the proximity of 
the charged particles in C-ST-1 arranged along the polymers in the static microgel provided 
enough counter ion repulsion to lessen the avaialble space for counter ion motion and thus 
decrease the overall counter ion percentage in the sphere going from the C-LS to the C-ST 
calculations. 

Although excluded volume had no effect on the configuration of the lattice sphere sim- 
ulations as described in section IIII B\ the specific arrangement of the excluded volume, in 
this case, arranged in a polymer chain, can slightly lessen the available space for counter ion 
condensation and thus lower the fraction of counter ions inside the microgel. 

V. MICROGEL IN A SPHERICAL CELL (C-MG) AND IN A PERIODIC BOX 
(P-MG) 

A. C-MG and P-MG Methodology 

After removing assumptions [Hand [2] to the lattice sphere (C-LS) and assumption [3] to the 
static micro gel (C-ST), we then removed the next assumption IH which looks at the effect 
of thermal fluctuations of the charged particles in the polymer. This is basically a microgel 
simulation, where the particles of the microgel are allowed to move, inside a spherical cell. 
Only the central particle of the microgel was kept fixed so as to keep the microgel in the 
center of the sphere. The parameters are basically the same as those of the microgel, and 
are listed in Table [TTl Aside from thermal fluctuations, the other difference of between this 
system and static microgel in a sphere (section IIV|) is the decrease in microgel radius, thus 
decreasing the free volume in between the chains. The parameters are listed in Table [III 
This simulation is labeled C-MG. As with the C-LS and C-ST simulations, the C-MG 
simulation was also done in a spherical cell model, with the electrostatics calculated as a 
direct sum. 

To determine the effect of the last assumtion, that of periodic boundary conditions and 
differences in electrostatic algorithm, we compared the results of C-MG to P-MG. The 
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methodology and parameters for the P-MG simulation has already been discussed in section 

Due to the thermal fluctuation of the entire microgel for both C-MG and P-MG, the 
radius of the microgel was also fluctuating, but generally maintaing its overall spherical 
shape. In order to determine the average radius of the microgel, we used the following 
equation^^ that relates the square radius of gyration of a rigid sphere to the radius of the 
sphere 

Rl = (10) 

where Rg is square radius of gyration of the sphere with radius r. 

B. Discussion: Thermal Fluctuations and Periodic Boundary Conditions 

The parameters and results (second to the last line) of the C-MG and P-MG are listed 
in Table IIIII The integrated counter ion probability is shown in figure [3] and the counter ion 
density is shown in figure HI Allowing the microgel particles to move allowed the microgel 
shrink from tq = 18 for the static microgel simulation (C-ST-1) to an average of tq = 14.96 
for C-MG and tq = 14.99 for P-MG. However, the fraction of counter ions within these 
radii for each simulation increased to 54.03% for C-MG and 53.95% for P-MG. Due to 
the difference of radius between PB-1 and these two microgel simulations, another PB 
calculation, PB-2 was done with a radius of ro = 14.96. The corresponding fraction of 
counter ions in the inner sphere for PB-2 also increased to 55.59%, resulting in a decrease 
of 2.6% when going from PB calculation to MG simulation of the same radius tq. As already 
seen in the discussion of the LS results in section IIIIBt the simluation results correspond to 
the PB calculation, with a slight decrease of a similar value (3.0%) in the fraction of counter 
ions in the lattice sphere. This shows that the thermal fluctuations of charged particles in 
a polymer do not affect the counter ion fraction inside the MG. 

The decrease of the radius of the charged microgel/sphere leads to an increase in fraction 
of counter ions in the micogel/sphere. This can be explained as follows. The decrease of the 
microgel/sphere radius, while keeping the total charge constant, corresponds to an increase 
in the overall charge density within the sphere, thus attracting more counter ions into the 
sphere. Our results of the PB calculations show that the gain in energy outweighs the loss 
of entropy, which the confined confined counter ions experience. 
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Both C-MG and P-MG have practically the same results — configuration as measured 
by Rg, fraction of counter ion inside the microgel, integrated counter ion probability (firgure 
E]) and counter ion densities (figure Hj). Thus the spherical cell model simulation, whose 
electrostatics is calculated as a direct sum, is equivalent to the corresponding simulation 
under periodic conditions, as long as the sytsems are set up correctly. This means that 
the system size is done ensuring enough electrostatic screening (as calculated by the Debye 
length) in order to satisfy the "isolated system" assumption of the cell model, and that the 
P3M electrostatics is tuned correctly. 

VI. CONCLUSIONS 

We have simulated a model polyelectrolyte microgel setup with the corresponding counter 
ions in a periodic box, using the results mainly to calculate the fraction of counter ions inside 
the microgel. We compared these results to a Poisson-Boltzmann numerical solution2^>2i in 
order to see how the mean field solution differs from the results of the simulation. We 
outlined the differences between these two models, listed as assumptions in section III CI 
and proceeded to perform a series of simulations aimed at removing these assumptions in a 
step-wise fashion in order to isolate the effects of each assumption. 

The very first assumption of going from a charged sphere with uniform distribution in the 
PB solver to having discrete charges showed a slight decrease in the fraction of counter ions 
inside the microgel. The fact that the sucessive "smearing out" of the discrete charges, thus 
approacing the uniform charge distribution in the PB solver, resluted in a sucessive approach 
of the LS results to that of the PB calculation shows that the discretization of charges causes 
the slight decrease of the fraction of counter ions inside the microgel. A further decrease is 
shown when going from the lattice sphere with evenly spaced charge distribution within the 
sphere to the static microgel configuration, where the charges are not as evenly distributed 
since they were lined up within the polymer. Aside from the discretization of charges, the 
density of charges of the microgel particles also plays a role in fraction of counter ions in the 
microgel, the two being directly proportional. 

Even with the many assumptions in between the MG simulations and PB theory, we 
have seen that the results of the two methods are quite close to one another. We can 
therefore safely assume that if PB models with configurations that follow the more complex 
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configurations found in molecular simulations of real systems are constructed, the results in 
terms of counter ion condensation from both approaches would be comparatively close, as 
long as we stick the the basic assumptions used in this study such as movable charges and 
counter ions, good solvent, and relatively isolated systems. This allows one to use the PB 
cell model to predict structural or thermodynamic properties of microgels in solution. 
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TABLE I: Parameters common to all simulations in this paper. These are: Lattice Sphere(C-LS), 
Static Microgel (C-ST), Microgel in a Spherical Cell (C-MG), and Microgel in a Periodic Box 
(P-MG). 



Parameter value 

ksT 1.0 

eij 1.0 

Ib 1.0 

unit mass m 1.0 

time step 0.012 

ensemble NVT 

thermostat Langevin 

friction coefficient (F) 1.0 

time steps for every snapshot 1000 
minimum number of snapshots for results 15000 

program used ESPResSo 
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TABLE II: Simulation parameters for the calculations discussed in sections IIIIBI and IIVBI The 
following parameters are common for all these calculations: the radius of the inner sphere tq = IScr 
whose total charge is qsp = +250.0, the total number of counter ions is nci = 250, the charge 
for each coutner ion is qci = —1.0, each counter ion is treated as a point particle with no aij 
interaction between them, all use the cell model and the electrostatic interactions are calculated 
via direct sum, the outer radius of the cell is i? = 52(T. The entries left blank in this table are 
parameters that are not applicable to the specific computation. 



Simulation label 


PB-1 C-LS-1 C-LS-2 C-LS-3 


C-LS-4 C-LS-5 


C-ST-1 C-ST-2 


number of charged particles {CP) 


250 


484 


894 


484 


894 


250 


250 


number of uncharged particles 

















9 


239 


excluded volume (ctlj) of CP 


1.0 


0.516 


0.279 


0.802 


0.654 


1.0 


1.0 


charge of particle {qcp) 


1.0 


0.516 


0.279 


0.516 


0.279 


1.0 


1.0 


percentage C/ in microgel/sphere 


52.83 49.80 


53.85 


52.78 


53.36 


52.21 


46.28 


45.67 


error of percentage CI 


2.01 


1.99 


2.00 


2.00 


1.98 


2.12 


2.09 
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TABLE III: Simulation parameters for the calculations discussed in section IV B[ The following 
parameters are common for all these calculations: the total charge of each particle in the inner 
sphere is qsp = +250.0, the total number of counter ions is nci = 250, the charge of each counter ion 
is qci = —1.0, the electrostatics of C-MG was calculated via direct sum, while the electrostatics 
of P-MG calculation which uses periodic boundary conditions was calculated using P3M. The 
entries left blank in this table are parameters and results that are not applicable to the specific 
computation. 



Simulation label 


PB-2 C-MG P-MG 


radius of sphere (constant) 


14.96 






boundary (Cell or Periodic box) 


cell 


cell 


box 


box length L (in a) 






84 


outer radius R (in a) 


52 


52 




number of charged particles (CP) 




250 


250 


number of uncharged particles 




239 


239 


excluded volume {(Tlj) of CP 




1.0 


1.0 


charge of particle (qcp) 




1.0 


1.0 


excluded volume of {(tlj) of CI 






1.0 


average squared radius of gyration 




134.27 


134.76 


error in R^ 




4.62 


4.51 


radius of microgel (eqn. [TO]) 




14.96 


14.99 


percentage CI in microgel/sphere 


55.59 


54.03 


53.95 


error in percentage CI 




2.36 


2.32 
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Figure [H Snapshot of a microgel P-MG simulation (shown in perspective) in a periodic 
box. Figure a. shows how the starting structure was constructed, while figure b. shows 
one equilibrated structure. The microgel is composed of 46 polymers, each composed of 
10 monomers (represented as grey shaded beads). The polymer also contains an additional 
29 beads acting as crosslinks (represented as dark beads, larger radius is done only for 
emphasis), which connect the ends of either three or four polymers. Thus the microgel has 
a total of 489 beads. Half of the monomers and 20 crosslinks, totalling 250 beads, have a 
+1.0 charge, giving a total charge of +250.0. The box contains 250 counter ions (not all 
shown), each with a charge of -1.0, making the box have an overall neutral charge. 

Figure [21 Models discussed in section III CI to go from the Poisson-Boltzmann spherical 
cell model (PB) to the microgel simulation (P-MG). The dark objects are the particles in 
the lattice (b) or the monomers in the microgel (c-e), while the lighter colored objects are 
the counterions that are free to move within the boundaries. The numbers on the arrows 
indicate the assumptions taken when going from one model to the next, as discussed in 
section Hi CI a. Poisson-Boltzann spherical cell model (PB). b. Lattice sphere (C-LS). c. 
Static microgel (C-ST). d. Cell microgel (C-MG). e. Microgel (P-MG). 

Figure [31 Integrated ion probability of simulations listed in Tables [IT] and IIIIl also corre- 
sponding to the counter ion densities in figure [H The inset is a magnification of the regions 
around tq = 18 (radius corresponding to table [TTl) and around tq = 14.99 (radius corre- 
sponding to table imp . The dots in the inset indicate the numbers reported at the last line 
in Tables [Hand [ml 

Figure [4l Counter ion densities of the simulations listed in Tables [IT] and IIII| also corre- 
sponding to the integrated counter ion probabilities in figure [31 The counter ion density 
graph of C-LS-4 is very simlar to C-LS-2, that of C-LS-5 to C-LS-3, that of C-ST-2 to 
C-ST-1, and that of C-MG to P-MG. These are therefore not included in this figure. 

Figure \5\ Arrangement of Excluded Volume. Solid circles represent the charged particles, 
open circles the counter ions, the grey shaded circle represents an uncharged particle, and 
the dashed circles the excluded volume cxij between the particles and the counter ions. a. 
Lattice sphere C-LS-1. The charged particles are at a distance of 4.81 apart for the C-LS-1 
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case, 3.67 for the C-LS-2 case, and 3.04 for the C-LS-3 case. After subtracting the two 
excluded volume radii from the charged particle, the distance x is still greater than Bjerrum 
length Ib-i thus allowing the counter ion to move freely around the dashed circle, b. Static 
microgel C-ST-1 without uncharged particles. Since the distance between charged particles 
is only around 2.0, there is no extra space for the counter ions to freely move once the 
excluded volume has been taken into account, c. Static microgel C-ST-2 with uncharged 
particles in between. The decrease in free space in b and c results in the lowering of the 
counter ion fraction for the ST case as compared to the LS. 
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